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Abstract 

A non-linear shell theory that includes transverse shear strains and its implementation in the 
material point method framework are discussed. The applicability of the shell implementation to model 
large deformations of thin shells is explored. Results suggest that an implicit time stepping scheme may 
be required for improved modeling of thin shells by the material point method. 

1 Shell Formulation 

The continuum-based approach to shell theory has been chosen because of the relative ease of 
implementation of constitutive models in this approach compared to exact geometrical descriptions of the 
shell. In order to include transverse shear strains in the shell, a modified Reissner-Mindlin assumption is 
used. The major assumptions of the shell formulation are d|2]| 

1. The normal to the mid-surface of the shell remains straight but not necessarily normal. The direction 
of the initial normal is called the "fiber" direction and it is the evolution of the fiber that is tracked. 

2. The stress normal to the mid-surface vanishes (plane stress) 

3. The momentum due to the extension of the fiber and the momentum balance in the direction of the 
fiber are neglected. 

4. The curvature of the shell at a material point is neglected. 

The shell formulation is based on a plate formulation by Lewis et al. 0). A discussion of the formulation 
follows. 

The velocity field in the shell is given by 

w(a, j3) = u(a, f3) + z uj{a, (3) x n(a, [3) + z n(Q!, [3) (1) 

where is the velocity of a point in the shell, u is the velocity of the center of mass of the shell, n is the 
normal or director vector, u> is the angular velocity of the director, (a, j3) are orthogonal co-ordinates on 
the mid-surface of the shell, z is the perpendicular distance from the mid-surface of the shell, and i is the 
rate of change of the length of the shell director. 
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Since momentum balance is not enforced for the motion in the direction of the director n, the terms 
involving z are dropped in constructing the equations of motion. These terms are also omitted in the 
deformation gradient calculation. However, the thickness change in the shell is not neglected in the 
computation of internal forces and moments. Equation Q can therefore be written as 

w(a, /3) = u(q, /?) + z r(a, /5) (2) 

where r, the rotation rate of n, is a vector that is perpendicular to n. 

The velocity gradient tensor for w is used to compute the stresses in the shell. If the curvature of the shell 
is neglected, i.e., the shell is piecewise plane, the velocity gradient tensor for w can be written as 



Vw 



+ r O n (3) 



where r n represents the dyadic product, and V'^'*) is the in-surface gradient operator, defined as, 

V(^) = [V( )] . l(^) . (4) 

The • represents a tensor inner product and I^*) is the in-surface identity tensor (or the projection operator), 
defined as, 

l(^) = I - n O n. (5) 

It should be noted that, for accuracy, the vector n should not deviate significantly from the actual normal to 
the surface (i.e., the transverse shear strains should be small). 

The determination of the shell velocity tensor Vw requires the determination of the center of mass velocity 
u of the shell. This quantity is determined using the balance of linear momentum in the shell. The local 
three-dimensional equation of motion for the shell is, in the absence of body forces, 

V • cr = p a (6) 

where sigma is the stress tensor, p is the density of the shell material, and a is the acceleration of the shell. 
The two-dimensional form of the linear momentum balance equation Q with respect to the surface of the 
shell is given by 

V(") • (cr) = p a . (7) 

The acceleration of the material points in the shell are now due to the in-surface divergence of the average 
stress (cr) in the shell, given by 

(8) 



1 r 

h J h- '^^^^ 



where /i+ is the "thickness" of the shell (along the director) from the center of mass to the "top" of the 
shell, is the thickness from the center of mass to the "bottom" of the shell, and h = + . The point 
of departure from the formulation of Lewis et al. (21 is that instead of separate linear momentum balance 
laws for shell and non-shell materials, a single global momentum balance is used and the "plane stress" 
condition cTzz = is enforced in the shell stress update, where the subscript zz represents the direction of 
the shell director. 

The shell director n and its rotation rate r also need to be known before the shell velocity gradient tensor 
Vw can be determined. These quantities are determined using an equation for the conservation of angular 
momentum 1 3 1, given by 

V(') • M - n • (cr) • I^") = ^ p r (9) 
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where r is the rotational acceleration of n, p is the density of the shell material, and M is the average 
moment, defined as 



M: =l("). 



h+ 

(t{z) z dz 

-h- 



• . (10) 



The center-of-mass velocity u, the director n and its rate of rotation r provide a means to obtain the 
velocity of material points on the shell. The shell is divided into a number of layers with discrete values of 
z and the layer-wise gradient of the shell velocity is used to compute the stress and deformation in each 
layer of the shell. 

2 Shell Implementation for the Material Point Method 

The shell description given in the previous section has been implemented such that the standard steps of the 
material point method f?] remain the same for all materials. Some additional steps are performed for shell 
materials. These steps are encapsulated within the shell constitutive model. 

The steps involved for each time increment Ai are discussed below. The superscript n represents the value 
of the state variables at time n At while the superscript n + 1 represents the value at time (n + 1) At. Note 
that At need not necessarily be constant. In the following, the subscript p is used to index material point 
variables while the subscript v is used to index grid vertex variables. The notation J2p denotes summation 
over material points and J2v denotes summation over grid vertices. Zeroth order interpolation functions 
associated with each material point are denoted by Sp^l while first order interpolation functions are denoted 
by^S- 

1 . Interpolate state data from material points to the grid. 

The state variables are interpolated from the material points to the grid vertices using the contiguous 
generalized interpolation material point (GIMP) method |5|. In the GIMP method material points are 
defined by particle characteristic functions Xp(x) which are required to be a partition of unity, 

^Xp(x) = 1VxgO (11) 
p 

where x is the position of a point in the body $7. A continuous representation of the property /(x) is 
given by 

/(x) = Y,fp Xp(x) (12) 

p 

where fp is the value at a material point. Similarly, a continuous representation of the grid data is 
given by 

= ^c/i> (13) 

where 

Sy{x) = 1 V X G n . (14) 

ti 

To interpolate particle data to the grid, the interpolation (or weighting functions) Sp]l are used, 
which are defined as 

^p JUpiin 
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where Vp is the volume associated with a material point, fip is the region of non-zero support for the 
material point, and 

^412 = 1 V XpSf^p. (16) 

V 

The state variables that are interpolated to the grid in this step are the mass (m), momentum (mu), 
volume (V), external forces (f^"'), temperature (T), and specific volume (v) using relations of the 
form 

= '^mp S^]l . (17) 
p 

In our computations, bilinear hat functions were used that lead to interpolation functions Sp]l 
with non-zero support in adjacent grid cells and in the next nearest neighbor grid cells. Details of 
these functions can be found in reference ||5l. 

For shell materials, an additional step is required to inhabit the grid vertices with the interpolated 
normal rotation rate from the particles. However, instead of interpolating the angular momentum, the 
quantity pp = mpTp is interpolated to the grid using the relation 

Pv = Y.PpS^]l. (18) 

p 

At the grid, the rotation rate is recovered using 

r„ = Pv/m^ (19) 

This approximation is required because the moment of inertia contains h"^ terms which can be very 
small for thin shells. Floating point errors are magnified when nip is multiplied by /i^. In addition, it 
is not desirable to interpolate the plate thickness to the grid. 

2. Compute heat and momentum exchange due to contact. 

In this step, any heat and momentum exchange between bodies inside the computational domain is 
performed through the grid. Details of contact algorithms used my the material point method can be 
found in references El|6l|71. 

3. Compute the stress tensor. 

The stress tensor computation follows the procedure for hyperelastic materials cited in reference HI. 
However, some extra steps are required for shell materials. The stress update is performed using a 
forward Euler explicit time stepping procedure. The velocity gradient Vw at a material point is 
required for the stress update. This quantity is determined using equation Q. The velocity gradient 
of the center of mass of the shell ( Vu) is computed from the grid velocities using gradient weighting 
functions of the form 

^p Jf^pfin 

so that 

Vnp = ^u,VS^]l. (21) 
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The gradient of the rotation rate (Vr) is also interpolated to the particles using the same procedure, 
i.e., 

Vrp = J]r,V5W. (22) 



The next step is to calculate the in-surface gradients V'-'^Vp and V^^^rp. These are calculated as 



V^up 



Vup • (I 
: Vrp . (I 



(23) 
(24) 



The superscript n represents the values at the end of the n-th time step. The shell is now divided into 
a number of layers with different values of z (these can be considered to be equivalent to Gauss 
points to be used in the integration over z). The number of layers depends on the requirements of the 
problem. Three layers are used to obtain the results that follow. The velocity gradient Vwp is 
calculated for each of the layers using equation For a shell with three layers (top, center and 
bottom), the velocity gradients are given by 



Vw7 



V^'^Up + h+ V(^)rp 



p 

,bot 



p 

V(")u„ - h- V(")r, 



I r. 



The increment of deformation gradient (AF) in each layer is computed using 

AFp = At Vwp + I 
The total deformation gradient (F) in each layer is updated using 

F;+' = AFp . F^ 

— n+l 

where F^ is the intermediate updated deformation gradient prior to application of the "plane 
stress" condition. 

The stress in the shell is computed using a stored energy function (W) of the form 



(25) 
(26) 
(27) 

(28) 

(29) 



W 



1 



K 



1 



1) - In J 



1 



G [tr(b) 



(30) 



where K is the bulk modulus, G is the shear modulus, J is the Jacobian (J = det F), and b is the 
volume preserving part of the left Cauchy-Green strain tensor, defined as 



J-3F«F^ 



(31) 



The Cauchy stress then has the form 

1 



K J 



i)i+^ 

J J 



-tr(b) 



(32) 



The "plane stress" condition in the thickness direction of the shell is applied at this stage using an 
iterative Newton method. To apply this condition, the deformation gradient tensor has to be rotated 
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such that its (33) component is aligned with the (zz) direction of the shell. The rotation tensor is the 
one required to rotate the vector = (0, 0, 1) to the direction about the vector x n^. If 6 is 
the angle of rotation and a is the unit vector along axis of rotation, the rotation tensor is given by 
(using the derivative of the Euler-Rodrigues formula) 



where 



The rotated deformation gradient in each layer is given by 



R = cos 61 (I 


— a (> 


D a) + a 


(g)a 




" 


-as 




A = 


as 





-ai 




-a2 


ai 






sin 6* A 



TTirot 



R . f!+' . R^ 



(33) 



(34) 



(35) 



The updated stress (<t™ ) is calculated in this rotated coordinate system using equation (I32l i. Thus, 



(36) 



An iterative Newton method is used to determine the deformation gradient component F33 for which 

o 

the stress component 1733 is zero. The "plane stress" deformation gradient is denoted F and the stress 
o 

is denoted a. 

At this stage, the updated thickness of the shell at a material point is calculated from the relations 

»1 o 



h. 



n+l 







Fzz{+z) dz 







h. 



n.+l 



Fzz{-z) dz 



(37) 
(38) 



where and /iq ai^e the initial values, and ^"^^ ^n+i updated values, of and h , 

respectively. 

In the next step, the deformation gradient and stress values for all the layers at each material point 
are rotated back to the original coordinate system. The updated Cauchy stress and deformation 
gradient are 



T7in+1 



R^ .F.R 
R^ • • R 



(39) 
(40) 



The deformed volume of the shell is approximated using the Jacobian of the deformation gradient at 
the center of mass of the shell 



n+l 



taO jn+1 
Vp Jp 



(41) 



4. Compute the internal force and moment. 
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The internal force for general materials is computed at the grid using the relation 



P 

For shell materials, this relation takes the form 



yn+l 
P 



(42) 



n+l 



(43) 



In addition to internal forces, the formulation for shell materials requires the computation of internal 
moments in order to solve for the rotational acceleration in the rotational inertia equation To 
obtain the discretized form of equation @, the equation is integrated over the volume of the shell 
leading to L2J 

- 5^ [(m, . V^i) . + (n, . (a,) . 5(5] V,= (^Y1 ^P ) • (44) 

p \ p / 

The average stress over the thickness of the shell is calculated using equation ^ and the average 
moment is calculated using equation dTOl . The trapezoidal rule is used in both cases. Thus, 



'-n+l J-h, 



'n+l 



(z) dz 



n + l 



J-h 



r^+\z)zdz 



where 



l(^) = I - 



(45) 
(46) 

(47) 



These are required in the balance of rotational inertia that is used to compute the updated rotation 
rate and the updated director vector. The internal moment for the shell material points can therefore 
the calculated using 



n+l 



(48) 



In practice, only the first term of equation ( R^ is interpolated to the grid and back to the particles. 
The equation of motion for rotational inertia is solved on the particles. 

5. Solve the equations of motion. 

The equations of motion for linear momentum are solved on the grid so that the acceleration at the 
grid vertices can be determined. The relation that is used is 



u, 



1 



(49) 



where f^"' are external forces. 
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The angular momentum equations are solved on the particles after interpolating the term 

= (m;+' . V5« . (50) 

p 

back to the material points to get nip. The rotational acceleration is calculated using 

12 Vp 



nip hi 



Tp) . I(^) 



(51) 



6. Integrate the acceleration. 

The lineal" acceleration in integrated using a forward Euler rule on the grid, giving the updated 
velocity on the grid as 

= < + At (52) 

For the rotational acceleration, the same procedure is followed at each material point to obtain an 
intermediate increment 

AFp = At Yp (53) 

The factor nip hp in the denominator of the right hand side of equation (I51t makes the differential 
equation stiff. An accurate solution of the equation requires an implicit integration or extremely 
small time steps. Instead, an implicit correction is made to Arp by solving the equation [9J 

[I + /? (I - n^)] Arp = A?p (54) 



where Arp is the corrected value of Arp and 



6E /At^^ 



Vpnip \ h J 

which uses the Young's modulus E of the shell material. The intermediate rotation rate is updated 
using the corrected increment. Thus, 

Tp =r'^ + Arp . (56) 

7. Update the shell director and rotate the rotation rate At this stage, the shell director at each 
material point is updated. The incremental rotation tensor AR is calculated using equation d33b with 
rotation angle 9 = \r\At and axis of rotation 

n" X r 

- - P P (57) 



Irip X Tp 



n"+i = AR . n" . (58) 



The updated director is 

Ip I^J-K, lip 

In addition, the rate of rotation has to be rotated so that the direction is perpendicular to the director 
using, 

r^+i = AR.?;^\ (59) 
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8. Interpolate back to the material points and update the state variables. 

In the final step, the state variables at the grid are interpolated back to the material points using 
relations of the form 

u^+i = ^<+i4i) (60) 

V 

Steps 1 through 8 are repeated for the next time step. 

3 Results 

Three tests of the shell formulation have been performed on different shell geometries - a plane shell, a 
cylindrical shell, and a spherical shell. 

1. Punched Plane Shell 

This problem involves the indentation of a plane, circular shell into a rigid cylindrical die of radius 8 
cm. The shell is made of annealed copper with the properties and dimensions shown in Table [2 



Table 1: Circular plane shell properties and dimensions. 



Po 


K 


G 


Thickness 


Radius 


Velocity 


(kg/m3) 


(GPa) 


(GPa) 


(cm) 


(cm) 


(m/s) 


8930 


136.35 


45.45 


0.3 


8 


100 



Snapshots of the deformation of the shell are shown in Figure [2 Substantial deformation of the shell 
occurs before particles at the edges tend to tear off. The tearing off of particles is due to the presence 
of large rotation rates (r) which are due to the stiffness of the rotational acceleration equation dSTl . 
The implicit correction does not appear to be adequate beyond a certain point and a fully implicit 
shell formulation may be required for accurate simulation of extremely large deformations. 

Particles in the figure have been colored using the equivalent stress at the center-of-mass of the shell. 
The stress distribution in the shell is quite uniform, though some artifacts in the form of rings appear. 
An implicit formulation has been shown to remove such artifacts in the stress distribution in 
membranes 1 10|. Therefore, an implicit formulation may be useful for the shell formulation. Another 
possibility is that these artifacts may be due to membrane and shear locking, a known phenomenon 
in finite element formulations of shells based on a continuum approach jTl|^. Such locking effects 
can be reduced using an addition hour glass control step 1 1 1 in the simulation. 

2. Pinched Cylindrical Shell 

The pinched cylindrical shell is one of the benchmark problems proposed by MacNeal and 
Harder 1 12 1. The cylindrical shell that has been simulated in this work has dimensions similar to 
those used by Li et al. fT^ . The shell is pinched by contact with two small rigid solid cylinders 
placed diametrically opposite each other and located at the midpoint of the axis of the cylinder. Each 
of the solid cylinders is 0.25 cm in radius, 0.5 cm in length, and moves toward the center of the 
pinched shell in a radial direction at 10 ms~^. The material of the shell is annealed copper (properties 
are shown in TableQ. The cylindrical shell is 2.5 cm in radius, 5.0 cm long, and 0.05 cm thick. 
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Figure 1: Deformation of punched circular plane shell. 
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Figure 2: Deformation of pinched cylindrical shell. 
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Snapshots of the deformation of the pinched cylindrical shell are shown in Figure|2l The deformation 
of the shell proceeds uniformly for 60 ms. However, at this time the increments of rotation rate begin 
to increase rapidly at each time step, even though the velocity of the center-of-mass of the shell still 
remains stable. This effect can be attributed to the stiffness of the rotational inertia equation. The 
effect is that extremely large rotation rates are produced at 70 ms causing high velocities and eventual 
numerical fracture of the cylinder. The problem may be solved using an implicit shell formulation. 

3. Inflating Spherical Shell 

The inflating spherical shell problem is similar to that used to model lipid bilayers by Ayton et 
al. fT4ll . The shell is made of a soft rubbery material with a density of 10 kg m^^, a bulk modulus of 
60 KPa and a shear modulus of 30 KPa. The sphere has a radius of 0.5 m and is 1 cm thick. The 
spherical shell is pressurized by an initial internal pressure of 10 KPa. The pressure increases in 
proportion to the internal surface area as the sphere inflates. 

The deformation of the shell with time is shown in Figure |3l The particles in the figure are colored 
on the basis of the equivalent stress. Though there is some difference between the values at different 
latitudes in the sphere, the equivalent stress is quite uniform in the shell. The variation can be 
reduced using the implicit material point method ITSll . 

4 Summary 

A shell formulation has been developed and implemented for the explicit time stepping material point 
method based on the work of Lewis et al. Q. Three different shell geometries and loading conditions have 
been tested. The results indicate that the stiff nature of the equation for rotational inertia may require the 
use of an implicit time stepping scheme for shell materials. 
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